program principal
interface
function leerFichero(A,vector)result(n) 
	implicit none
	integer::n,j,i
	real(kind=8),dimension(:,:),allocatable::A
	real(kind=8),dimension(:),allocatable::vector
end function leerFichero

subroutine substiDelante(V,b,dimen)
	real(kind=8),dimension(:,:),allocatable::V
	real(kind=8),dimension(:),allocatable::b
	integer::dimen
	integer::k
end subroutine

subroutine escribe(LU,dimen,x)
	real(kind=8),dimension(:,:),allocatable::LU
	real(kind=8),dimension(:),allocatable::x
	integer::dimen,k
end subroutine

subroutine substiDetras(V,result,dimen)
	real(kind=8),dimension(:,:),allocatable::V
	real(kind=8),dimension(:),allocatable::result
	integer::dimen
	integer::k
end subroutine

subroutine lu(dimen,V)
	real(kind=8),dimension(:,:),allocatable::V
	integer::i,j
	integer::dimen
end subroutine lu

subroutine jacobi(A,b,dimen,it,tol)
	real(kind=8),dimension(:,:),allocatable::A
	real(kind=8),dimension(:),allocatable::b,z,x,r
	real(kind=8)::norma2r,norma2b,tol
	integer::dimen,it,i,k
end subroutine

subroutine relajacion(A,b,dimen,it,tol)
	real(kind=8),dimension(:,:),allocatable::A
	real(kind=8),dimension(:),allocatable::b,z,x,r
	real(kind=8)::norma2r,norma2b,tol,w
	integer::dimen,it,i,k,j
end subroutine

subroutine ampliar(S,V,dimen,n)
	real(kind=8),dimension(:,:),allocatable::S,V,I
	integer::dimen,j,n
end subroutine

subroutine datos(it,tol)
		integer::it
		real(kind=8)::tol
end subroutine

function menu()
	integer::opcion
end function menu

end interface
integer::dimen,op,n
real(kind=8),dimension(:,:),allocatable::A,V
real(kind=8),dimension(:),allocatable::vector
integer::it
real(kind=8)::tol
	n= leerFichero(A,vector)
	dimen=n*3
	call ampliar(A,V,dimen,n)
	op= menu()

	select case (op)
			case (1) 
				call lu(dimen,V)
				call substiDelante(V,vector,dimen)
				call substiDetras(V,vector,dimen)
				call escribe(V,dimen,vector)
			case (2) 
				call datos(it,tol)
				call jacobi(V,vector,dimen,it,tol)
			case (3) 
				call datos(it,tol)
				call relajacion(V,vector,dimen,it,tol)
		end select


end program
	function leerFichero(A,vector) result(n)
		implicit none
		integer::n,j,i
		real(kind=8),dimension(:,:),allocatable::A
		real(kind=8),dimension(:),allocatable::vector
		open(unit=1, file='datos.dat',status='unknown')
		read(1,*)n
		allocate(A(n,n))
		allocate(vector(n*3))
		read(1,*)A
		read(1,*)vector
		A=transpose(A)
		close(1)
	end function
	
	subroutine substiDelante(V,b,dimen)
		real(kind=8),dimension(:,:),allocatable::V
		real(kind=8),dimension(:),allocatable::b
		integer::dimen
		integer::k
		do k=2,dimen,1
			b(k)=b(k)-dot_product(V(k,1:k-1),b(1:k-1))
		end do
	end subroutine

	subroutine escribe(LU,dimen,x)
		real(kind=8),dimension(:,:),allocatable::LU
		real(kind=8),dimension(:),allocatable::x
		integer::dimen,k
		do k=1,dimen,1
			print*,x(k)
		end do
		open(unit=2, file='LU.dat',status='unknown')
		do k=1,dimen,1
			write(2,*)LU(k,1:dimen)
		end do
		close(2)
	end subroutine

	subroutine substiDetras(V,result,dimen)
		real(kind=8),dimension(:,:),allocatable::V
		real(kind=8),dimension(:),allocatable::result
		integer::dimen
		integer::k
		result(dimen)=result(dimen)/V(dimen,dimen)
		do k=dimen-1,1,-1
			result(k)=(result(k)-dot_product(V(k,k+1:dimen),result(k+1:dimen)))/V(k,k)
		end do
		
	end subroutine



	subroutine lu(dimen,V)
		real(kind=8),dimension(:,:),allocatable::V
		integer::i,j
		integer::dimen
		do i=1,dimen,1
			V(i,i)=V(i,i)-dot_product(V(i,1:i-1),V(1:i-1,i))
			do j=i+1,dimen,1
				V(i,j)=V(i,j)-dot_product(V(i,1:i-1),V(1:i-1,j))
				V(j,i)=(V(j,i)-dot_product(V(j,1:i-1),V(1:i-1,i)))/V(i,i)
			end do
		end do
	end subroutine

	subroutine jacobi(A,b,dimen,it,tol)
		real(kind=8),dimension(:,:),allocatable::A
		real(kind=8),dimension(:),allocatable::b,z,x,r
		real(kind=8)::norma2r,norma2b,tol
		integer::dimen,it,i,k
		
		allocate(x(dimen))
		allocate(z(dimen))
		allocate(r(dimen))

		open(unit=3, file='Jacobi.dat',status='unknown')

		x(1:dimen)=0

		write(3,*)'iteraciones: ',it,'tolerancia: ',tol
		write(3,*)'vector inicial: '
		write(3,*)x(1:dimen)
		write(3,*)
		do k=1,it,1
			z=MATMUL(A,x) 
			r=b-z
	
			norma2r=sqrt(dot_product(r,r))
			norma2b=sqrt(dot_product(b,b))

			if(norma2r<tol*norma2b) then 
				exit
			end if

			do i=1,dimen,1
				z(i)=r(i)/A(i,i)
				x(i)=x(i)+z(i)
			end do

			write(3,*)'iteracion: ',k
			write(3,*)'aproximacion: ',x(1:dimen)
			write(3,*)'error relativo: ',norma2r/norma2b
			write(3,*)
		end do	
		close(3)
	end subroutine

	subroutine relajacion(A,b,dimen,it,tol)
		real(kind=8),dimension(:,:),allocatable::A
		real(kind=8),dimension(:),allocatable::b,z,x,r
		real(kind=8)::norma2r,norma2b,tol,w
		integer::dimen,it,i,k,j

		allocate(x(dimen))
		allocate(z(dimen))
		allocate(r(dimen))
		print*,'¿que w desea?: '
		read(5,*)w

		open(unit=4, file='Relajacion.dat',status='unknown')
		x(1:dimen)=0

		write(4,*)'iteraciones: ',it,'tolerancia: ',tol
		write(4,*)'vector inicial: '
		write(4,*)x(1:dimen)
		write(4,*)

		do k=1,it,1
			z=MATMUL(A,x) 
			r=b-z

			norma2r=sqrt(dot_product(r,r))
			norma2b=sqrt(dot_product(b,b))

			if(norma2r<tol*norma2b) then 
				exit
			end if

			z(1)=r(1)/A(1,1)
			do i=2,dimen,1
				z(i)=(r(i)-w*(dot_product(A(i,1:i-1),z(1:i-1))))/A(i,i)
			end do

			x=x+w*z

			write(4,*)'iteracion: ',k
			write(4,*)'aproximacion: ',x(1:dimen)
			write(4,*)'error relativo: ',norma2r/norma2b
			write(4,*)
		end do
		close(4)
	end subroutine

	subroutine ampliar(S,V,dimen,n)
		real(kind=8),dimension(:,:),allocatable::S,V,I
		integer::dimen,j,n
		allocate(V(dimen,dimen))
		allocate(I(n,n))
		I=0
		do j=1,n,1
			I(j,j)=-1
		end do

		V=0
		V(1:3,1:3)=S
		V(1:3,4:6)=I
		V(4:6,1:3)=I
		V(4:6,4:6)=S
		V(4:6,7:9)=I
		V(7:9,4:6)=I
		V(7:9,7:9)=S

	end subroutine

	subroutine datos(it,tol)
		integer::it
		real(kind=8)::tol
		print*,'Cuantas iteraciones desea?'
		read(5,'(i2)')it
		print*,'Cual es la tolerancia de error?'
		read(5,*)tol
	end subroutine
	
	function menu()result(opcion)
	integer::opcion
		print*,'Que factorizacion desea?'
		print*,'1.- LU'
		print*,'2.- Jacobi'
		print*,'3.-Relajacion'
		read(5,'(i2)') opcion
		
	end function


	
